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Abstract 

We present new sets of next-to-leading order fragmentation functions describing the production 
of charged pions, kaons and protons from the gluon and from each of the quarks, obtained by fitting 
to all relevant data sets from e + e _ annihilation. The individual light quark flavour fragmentation 
functions are obtained phenomenologically for the first time by including in the data the light 
quark tagging probabilities obtained by the OPAL Collaboration. 
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I. INTRODUCTION 



Theoretical predictions for future experiments are necessary for determining the kine- 
matic regions of validity of the Standard Model (SM). Such predictions depend on constants 
which must be determined from past experiments since these quantities are otherwise uncal- 
culable, either because no theory exists which can determine them from more fundamental 
parameters, or because the solutions of the current theory are insufficient to determine them 
from the SM parameters. 

Quantum Chromodynamics (QCD), the theory of the strong interaction and one of the 
theories that make up the SM, is required in the description of processes involving hadrons. 
The best tool for solving QCD to perform such descriptions is perturbation theory How- 
ever, perturbative QCD (pQCD) can only describe the high energy components of the cross 
section, while a process will contain low energy components if a hadron is in the initial 
state or is observed in the final state. Fortunately, from the Factorization Theorem, the 
low and high energy scale components of such processes can be separated. The low energy 
components are universal and so can be used to make predictions. Since they cannot yet be 
reliably calculated from QCD, they must be extracted from experimental data. 

The pQCD description of data involving the inclusive production of hadrons requires 
fragmentation functions (FFs), which form the low energy components of such processes 
and describe the inclusive emission of a hadron from a quark or gluon (parton) for every 
momentum fraction. One reason FFs are important is that model independent predictions 
of LHC cross sections in which a hadron is detected in the final state depend on them. There 
are many theoretical obstacles to the extraction of FFs from data: The Dokshitzer-Gribov- 
Lipatov-Altarelli-Parisi (DGLAP) [1] evolution equation for FFs is only known to next-to- 
leading order (NLO), and is furthermore unreliable at small and possibly even intermediate 
momentum fractions of the emitted parton, where the only reliable determination of FFs is 

Q 

via the Modified Leading Logarithm Approximation (MLLA) |2j. Despite these problems, 
FFs at intermediate to large momentum fractions obtained from fits to data now yield 
compatible results with other data sets Q]. 

Much precise data from e + e~ colliders now exists for the production of the three lightest 
charged hadrons, which are the pion (7r ± ), kaon (i^ ± ) and proton (p/p). In much of this 
data, the observed hadron is identified as one of these particles, and the emitting parton 
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is identified as either a gluon, light (u, d and s) quark, c quark or b quark, which allowed 
for a precise determination of the corresponding individual FFs in Refs. P, Q| . How- 
ever, the individual light quark FFs could only be extracted by making reasonable physical 
assumptions. 

Since this analysis, the OPAL Collaboration has presented light flavour separated mea- 
surements on light charged hadron production p for the e + e~ centre-of-mass (CM) energy 
y/s = Mz, allowing for the first time the extraction of flavour dependent FFs of light quarks. 
In this OPAL analysis, high energy mesons (tt ± , K and Kg) and baryons (p/p and A) were 
identified in the large Z boson decay data sample and used as tagging products. In addi- 
tion, high momentum e^, /i ± and D*^ particles and identified bottom events were used to 
measure heavy flavour backgrounds in the above meson and baryon sample. As suggested in 
Ref. [jj and precisely studied in a recent analysis by the SLD Collaboration [8], these high 
energvparticles give information about the original quark. For more details see the OPAL 
work j6 : j, where it is explained how the Collaboration measured the probability T]^(x p , s) for 
a quark flavour a to develop into a jet containing the particle h with a momentum fraction 
x larger than x p = 2ph/\fs. 

Since the valence structure of the proton is uud, knowing the difference between the 
individual light flavour FFs, in particular for u and d quarks into K^, is very much needed 
for predicting the inclusive cross sections for the productions of these hadrons in collisions 
involving protons, such as ep, pp and pp collisions. For example, results from the inclusive 
production of hadrons in pp collisions provide the baseline to which one compares heavy- 
ion collision results in order to determine the properties of the hot quark-gluon plasma 
Tests presented in Ref. of the KKP FFs in the process p + p — > + X, where 
h ± are light charged hadrons, were generally successful, as was a recent check of the pion 
FFs by comparison to p + p — > n° + X data (taking 7r° = ~ (n + + 7T~)) from the PHENIX 
Collaboration at RHIC. However, it is likely that the inaccuracy on the information on 
the u, d and s quark FFs canceled out due to the superimposition of the hadrons in k* 1 



In this paper, we update the analysis of Ref. [3j by including the data of Ref. Jf| in 
the fit to obtain for the first time a phenomenological determination of the individual light 
quark FFs for each light charged hadron species. Since we do not impose those physical 
assumptions on the light quark FFs that were used in Ref. in our calculation of the cross 
sections used for the fit, the other FFs extracted in this fit are also more reliable. In Section 
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HTl we summarize the basic theoretical tools used in our calculations for the fit. In Section ITTT1 
we justify specific choices for our fit such as the data used and the FF parameterization. Our 
results are then presented in Section HVI and finally in Section we present our conclusions. 
The details of the longitudinal cross section calculation are given in Appendix [X] 



II. FORMALISM 



The optimal way to determine FFs is to fit them to measurements of the processes 
e + + — > (7, Z) —>■ a + a — > h + X , where a is the tagged quark, h is a detected hadron and 
X is the remaining unobserved part of the final state. In a typical experiment the hadron is 
only detected if its species h belongs to a specified set of hadron species Sh and the species 
of the tagged quark a belongs to a set of flavours S^. Writing the CM momentum of the 
observed hadron as X\fs/2, the data for such a process are typically presented as 

The total cross section a a is given to NLO by 

*a{s) = a (s)N c Q a (s) (1 + 2a s {s)) , (2) 

where cr = 47ra 2 /(3s) is the leading order (LO) cross section for the process e + + e~ — > 
7 — > n + + N c is the number of colours and a s (/i 2 ) = a s (fj) /(2tt). Q a (s) is the effective 
electroweak charge of quark a [l^ . 

From the Factorization Theorem, the higher twist component of the differential cross 
section in Eq. ((TJ is oi 0(1/^) or less and may and will be neglected in this paper, while the 
leading twist component is obtained by convoluting the corresponding high energy partonic 
cross sections with the FFs D%(y, Mjf), where y is the fraction of the momentum of parton 
a taken away by the produced hadron h and Mf is the factorization scale. This may be 
written concisely by taking y = x/z, in which case 



do^ , f l dz 



drr a 

aa a,NS 

dz 



(3) 
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where, for the emission of a hadron h from a quark a, the non-singlet partonic cross section 
da^s/dz contains only and all those contributions from diagrams in which the quark line 
connected to the electroweak vertex and the quark line emitting the hadron h are the same, 
while the pure singlet partonic cross section da% PS /dz contains all other contributions. Since 
the Z boson only splits into a quark a and its antiquark a, each partonic cross section is 
proportional to Q a , and thus may be written 46] 



dz 



(z, s, Mj) =a (s)Q a (s)C NS I z, a s (s), In 
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a,PS 



dz 
daft 
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PS 
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and 



(4) 



dz 



f (z, s, Mj) =2a (s)Q a (s)C g z, a s (s),\n -L 



Mi 



where the Cj are the perturbatively calculable coefficient functions. Uf is the number of 
active quark flavours. For the choice M| = s, the Ci(z, a s (s), 0) = Ci(z,a 8 (s)) for the 
unpolarized (i.e. summed over transverse and longitudinal components) cross section are 
given to NLO by Q 

" '2tt 2 

' ~ 2 
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(5) 

(6) 
(7) 



Note that the pure singlet contribution only enters at NNLO. In contrast, in the longitudinal 
cross section 



(8) 



there is a contribution from the pure singlet sector at NLO, while the gluon FF enters at 
LO (see Appendix IX]) . 

It is clear that we only apply electroweak theory to LO. We can therefore easily see 
that Eq. (JHJ) for the cross section when quark a is tagged is a physical observable, since 
it can be obtained by differentiating the untagged cross section (Eq. (jHJ) with a summed 
over all flavours) with respect to \nQ a , where Q a is the effective electroweak charge of 
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quark a discussed above. Therefore the tagged cross section is formally independent of 
the factorization and renormalization scales and schemes, as it must be to qualify as an 
observable. 

M 2 

For Mj ^ s, the coefficient functions will contain terms of the form a"(s) ln p where 
p — n, n — 1, which will spoil the convergence of the series unless Mj = O(s). Thus, in 
order to be able to describe data over a large range in s, the dependence of the FFs on Mj 
must be known. Fortunately this can be calculated using the DGLAP equation, 

j^pD^, Mj) = W' d -l Pab (i a,(Mj)) flf (v, Mj), (9) 

/ b=g,q Jz y \y / 

where the a — > b splitting functions P a b are perturbatively calculable, and are known to 
NLO. Therefore, in the calculation of cross sections it is sufficient to know the FFs at just 
one factorization scale Mf = M . 

The DGLAP equation is however not valid when z is small, since due to soft gluon 
emission the P a g(z, a s ) contain terms which behave in the limit z — > like (a™/z) i n 2n ~ 1 ~ m ^ 
where m — 1, ...,2n — 1 labels the class of terms (finite terms which behave like a™ when 
integrated over the range < z < 1 are classified as m = 2n), and are therefore unreliable 
in this limit. This implies that the cross section cannot be reliably calculated at small x, 
and the FFs (z, Mq ) cannot be fitted at small z. In this case a description of the data 
requires an alternative approximation such as the MLLA, which is beyond the scope of this 
paper. 

Dependence on the factorization scale is introduced in the usual way. Specifically, the 
FFs are evolved to Mj = kfS, where kf is a constant which is taken to be equal to 1 for 
the main fit, and 1/4 and 4 in two further fits to determine the theoretical errors on fitted 
parameters. We counter-balance this kf dependence at NLO using the result (where the x 
dependence, integrals, discrete labels, sums and charges have been removed for brevity) 

C(a s (s), In kf) = C(a,(s), 0) - In k f C(a s (s), 0)P(a s (s)). (10) 

Dependence on the renormalization scale \i is introduced by choosing fj, 2 = ks, where A; is a 
constant chosen to obey k = kf. At NLO, this amounts to replacing a s (s) in the coefficient 
functions with a s (ks). 

The fastest and most accurate way of calculating a cross section is in Mellin space, defined 
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by the transformation 



Fsl(n, s) = [ dxx^F^ix, s), (11) 
Jo 

since convolutions such as that in Eq. © become simple products. In particular, Eq. © 
becomes 

<! D h a (n, M}) = P ab (n, a s {Mj))D h b {n, Mf), (12) 



din M 2 f 

1 b=g,q 

which can be solved analytically order by order. The cross section in x space can then be 
obtained numerically via the inverse Mellin transform, 

F*»(x,s) = ^-j c dnx- n F^(n,s), (13) 

where C is a contour in Mellin space from Im(n) = — oo to Im(n) = oo, which passes to the 
right of all poles. 

Predictions for data averaged over an x-bin in the range x\ < x < Xh are calculated from 
the formula 

(F°«)(x h x h , s) = r dxF§;(x, s). (14) 

%h -El J x t 

This integral over x can be done analytically in Eq. ([13)1. 

i i r r l— n _ T l-n 

{ Fl mw) = __ j^n^—^F^M, (15 ) 

giving a further advantage for working in Mellin space that no extra numerical integration 
is required to obtain x-bin averaged cross sections. 

The light flavour separated data in Ref. [3] may be interpreted as the probability for a 
tagged quark flavour a to inclusively emit a hadron of type h with momentum greater than 
Xpy/s/2, in which case the corresponding theoretical result for such data may be calculated 
from the formula 

V h a (x p , s) = I dxF$(x, s) = (1 - x p )(F$)(x p , 1, s), (16) 

J Xp 

and we note that for this expression the tq{x p ,s) constrain the FFs at large momentum 
fraction even more than the Fg^(x,s). However, the experimental definition of the rfc is 
a little more subtle. For a given number iV a of e + e~ annihilation events in which a quark 
a is tagged, the number N a ^h of times that an event hemisphere, defined to be the two 
regions separated by the plane perpendicular to the thrust axis for each event, contains a 
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particle h with x > x p is determined. Therefore, at LO, where a and a are never in the 
same hemisphere, r]^(x p , s) is given by the integral over D%(x, s) in the range x p < x < 1, 
and this result is consistent with Eq. At NLO the quark a can emit a gluon which in 

turn emits the hadron h according to the gluon FF Dg (see Eq. fljJJ)). In the measurement of 
f]a{x p , s), processes in which the gluon is in the opposite hemisphere from the quark a that 
emitted it are excluded. However, such processes contribute to Eq. (JTBJ). Fortunately, such 
events in which the gluon is emitted with a large angle with respect to the quark a are very 
rare and should contribute very little both to Eq. (fTBj) and the measured 77^. 



III. METHOD 

In this Section we describe our method for obtaining FFs from data. As in Ref. Q], where 
a detailed discussion of all available data sets is given which will not be repeated here, we 
use identified hadron data with and without flavour separation from DELPHI and SLD 
|8j, and identified hadron data without flavour separation from ALEPH and TPC Q. 
In addition, we use identified hadron data with flavour separation from TPC |l7j], which was 
used in Ref. J4] but not in Ref. js]]. Furthermore, for the first time we also include the light 
: lavour separated measurements of quark tagging probabilities from the OPAL Collaboration 



However, we exclude unidentified hadron data since, although such data is accurate, it 
is typically contaminated with charged particles other than the 7r ± , and p/p. Such data 
was used in Ref. |3[, leading to consistent results. However, since in this analysis we aim 
for more reliable FFs, we use only hadron species separated measurements. We also exclude 
data for which x\ < 0.1, since the prediction for the cross section is unreliable in this region 
as a result of the logarithms from soft gluon emission mentioned in Section [H] After fitting, 
we then compare cross sections calculated from our FFs and a s [Mz) with the unidentified 



hadron data with flavour separation from TPC 

from ALEPH |l8,l3, DELPHI 14 1 and OPAL 20j, without flavour separation from SLD 



13], with and without flavour separation 



2l| and OPAL [22] and 



the unidentified hadron gluon-tagged three-jet data from ALEPH 
;he identified hadron tagging probabilities with heavy quark flavour separation from OPAL 
f|. The latter data set is not included in the fit since the heavy quark FFs are much better 
constrained by the larger quantity and quality of heavy quark-tagged data from DELPHI, 
SLD and TPC. In all data, correlation effects between data points are not yet known, and 
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therefore in the calculation of the covariance matrix for the \ 2 to be minimized, we fix the 
off-diagonal elements to zero, i.e. we assume that the data points are uncorrelated and that 
the error on each one is given by its statistical and systematic errors added in quadrature. 
Note that this common deficiency in published data limits the reliability of results obtained 
from analyzing them. 

All theoretical quantities are calculated to NLO in the MS scheme. For our main fit, we 
evolve the FFs from My = M to Mf = a/s, and vary the scales as described in Section |H] 
to determine the theoretical errors on a s (Mz). We take M = v2 GeV for a = u,d, s, g. As 
Mf is increased from Mq to a/s, the number of flavours used in the evolution of the FFs and 
the strong coupling is first set to rif = 3 and only the light quark and gluon FFs are non zero 
until Mf = m(i] c ) = 2.9788 GeV, where the charm FF is set equal to its initial distribution 
and included in the set of FFs to be evolved, and the number of flavours is taken to be nf = 4. 
The bottom FF is treated in the same way, being introduced when Mf = m(T) = 9.46037 
GeV. Both flavour thresholds are respectively twice the pole masses of these two heavy 
quarks, and therefore perturbative matching conditions are required at NLO. Rather than 
implementing this matching explicitly, we define our heavy quark FFs to be the complete 
ones, not just the intrinsic FFs, which means the matching term, dependent on the gluon 
FF, is absorbed into them. Our FFs are summed over hadrons which are of the same species 
but opposite charges, and averaged over quark and antiquark. We do not consider cross 
sections which depend on the difference between quark and antiquark FFs summed over any 
given set of emitted hadrons, although it must be noted that this difference is zero when this 
set contains a sum over charges, by charge conjugation invariance. Since we use accurate 
data at 29 GeV and 91.2 GeV, we are in a position to extract the parameter a s (Mz), the 
quantity which determines the running of a s (/i 2 ). We therefore free this parameter in our 
fit. The matching on a s (/i 2 ) is implemented by determining A^ from a s (Mz) } and then 
using it to determine A^L and A^ from the NLO relations given in Ref. 1231 (these were 
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checked using the results of Ref. 



24j). We choose the usual parameterization 



D^(x,M$) = Nx a (l -xf (17) 

for each of our FFs. In Mellin space, the FFs are then proportional to T(n + a)/T(n + 
a + (3 + 1) ~ l/(n + a) for n ~ —a, and this behaviour persists even after evolution and 
convolution with coefficient functions. For such behaviour, the numerical evaluation of Eq. 
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(|15|) is best performed with the integration variable < t < 1 and contour defined through 

3 1 2(1 -i) hit 

n = cH r H h — ( , (18) 

3 + ln^ I -xi ln^ ' y ' 

where the real constant c is chosen such that the contour lies to the right of all poles, since 
as t — > 1 the integrand in the integral over t becomes a finite constant, while as t — > 
the integrand vanishes like exp((— 2(1 — i) hit/ \a{l/xi)) lnx). As a result of the second and 
third term in Eq. ()18j) . the intersection of the contour with the lm{n = 0) line goes from 



n = 1 to n = oo as i goes from to 1. This approximately follows the saddle point |25( of 
the integrand, thus ensuring the contour is close to the contour of steepest descent, which 
gives the fastest convergence of the integral. 

In Ref. jjj, no data was used which could allow for the difference between the d and s 
FFs to be determined. (The FFs for the u can be determined since its electroweak charge is 
different to that of d and s.) The authors constrained this difference by imposing the valence 
quark structure at all momentum fractions and SU(3) invariance, giving the relations 

Df(x,M 2 )=Df (x,M 2 ), 

Df?(x, Ml) = D^(x, Ml) and (19) 
D^(x,Ml)=2Df (x,M 2 ). 



Such constraints can be implemented by fixing the parameters N, a and (3 of the FFs on 
the right hand sides to be equal to those of the FFs on the left hand side, with the exception 
that the parameter N of Du must be fixed to twice the value of that of D P J^ 47 1. With 



such conditions on the parameterization, a good fit to the data used was obtained. 



The first line in Eq. (p 

expected to be accurate |2 



also follows from SU(2) isospin invariance, and is therefore 
. Indeed, the approximate result = rfc implied by this 
relation is found to hold within 2% for x p > 0.2. However, the second line in Eq. (|T9~|) is 
expected to be strongly violated since the s quark has a significantly larger mass than the 
u quark. Already in 1977, Field and Feynman j^l assumed that due to the larger mass of 
s quarks, the s —>■ K + transition should happen more frequently than the u —>■ K + one 
because less energy is needed for the creation of a uu pair from the vacuum than for a 
ss pair. This is measured by the suppression factor 7 S of strange quarks, which is known 
from various strange/non-strange hadron production rates to be around 7 S ~ 0.3. (For a 
compilation, see Ref. 



.) The third line in Eq. (fT9*j) . assumed earlier also in Ref. [28], 
10 



can also be justified for x — > 1 by the valence ratios and dimensional counting powers [29]. 
Indeed, in the OPAL analysis of Ref. (|, the ratio rfj* '/rju^ is consistent with 0.5 for all 
x p > 0.2, but only inside the rather large errors. However, decays from heavier baryons such 
as A or A resonances might change this ratio. Furthermore, within the LUND string model 
[sol ] the actual value of the ratio rfj^ jnfj^ at large x p would be a direct measure of the size 
of the suppression of diquarks with spin 1 relative to those with spin 0, since Fermi-Dirac 
statistics requires a uu diquark to have angular momentum L — 1. In summary, all relations 
in Eq. (|19|h particularly the last two, may be violated to a possibly relevant degree, but 
in any case since we will use the data of Ref. jfj in our analysis, we shall not impose any 
relations between the light quark flavour FFs. 



IV. RESULTS 

In this Section we report the results obtained from the fit described in Section ITTT1 We 
obtain 

a s (M z ) = 0.1176j»[exp]lJ:SS[theo] = 0.1176i£g°g. (20) 

This is equivalent to the result A^ = 221 ± 74[exp]^ [theo] MeV. The experimental errors 
are obtained by varying a s (Mz), keeping all other parameters fixed, until Xdf increases by 
unity. The theoretical errors, determined using the method described in Section |H1 turn out 
to be negligible relative to the experimental ones, most likely because the x range of the 
data used is very limited. The second result in Eq. ([20)1. whose upper and lower errors are 
obtained by adding the upper and lower errors respectively of both sources in quadrature, is 
consistent with the KKP result [31] of a s (M z ) = O.mO^oora (which includes the theoretical 
error). In Tabled we show the values of the remaining, FF parameters obtained from the 
fit. Since N and (3 are highly correlated and the large x data generally has the largest errors, 
for some FFs these two parameters are large. However, over the range 0.1 < x < 1, all FFs 
are of similar order in magnitude. Also shown in Table ^ is the symmetrized propagated 
experimental error on each parameter. This quantity is the average of the two resulting errors 
obtained by varying the parameter, keeping the other parameters fixed, until \ 2 increases 
by 1 from its minimum value. The correlated errors between the parameters are expected 
to be of a similar order of magnitude to the purely statistical errors shown. Note that these 
results show no obvious consistency with Eq. (|19|) . With the inclusion of correlation effects 
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in the data, a deeper investigation into parameter errors would be worthy. 

We obtained Xdf = 1-15, indicating an overall good description of the data. The resulting 
Xp F values for the OPAL light quark tagging probabilities from Ref. (J are shown in Table 
121 The description of the data in which K or p/p is detected is excellent, except for the 
process d —>■ p/p. For this and the 7r ± data, which has the highest accuracy, a fit without 
the data points at x p = 0.2 results in all values of Xdf being around unity, although the 
resulting FFs from that fit are not considerably different to those from the main fit. This 
data, together with the corresponding theoretical curves calculated from our FF set (labeled 

fi n 

AKK), and with the curves from the sets of Ref. \$ (labeled KKP) and Ref. [4| (labeled 
Kretzer), are shown in Fig. ^ We see that for the s,d — > transitions, the corresponding 
AKK curves are in good agreement with the data while the KPP and Kretzer curves strongly 
disagree. The Kretzer set fails to lead to a decent description for the rf d data, but otherwise 
all ^ data is well described by all three sets. Our set and the KKP set lead to a good 
description of the p/p data (which were not used in the determination of the Kretzer set). 
Fig. El shows the heavy quark tagging probabilities, which were not used in the fit, together 
with the corresponding theoretical curves from the same FF sets as were used in Fig. [0 In 
Table H3 we list the corresponding Xbf values. Clearly these values are unacceptably high. 
In order to check that this was not a result of the inadequacy of our parameterization to 
allow for a description of both small x and large x data (since, as discussed around Eq. 
(I16j) . the OPAL quark tagging probabilities provide more constraints on the FFs at large 
x), we performed three new fits which included the heavy quark tagging probabilities, the 
first being otherwise similar to the main fit, the other two having the following differences: 
For the second fit, the quark FFs were modified by multiplying the right hand side of Eq. 
()17|) with (1 + jx), with 7 different for each quark FF and fixed to zero for the gluon FF, 
and each 7 was included in the set of free parameters to be fitted. In the third fit, all 
x\ < 0.2 data were excluded. No significant improvement to the description of the heavy 
quark tagging probabilities was obtained in all three fits. We therefore assume that this 
discrepancy is caused by the inclusion of large angle gluon emission effects in Eq. (116)1 . as 
described at the end of Section |Hj However, since we have sufficient data to constrain the 
heavy quark FFs, we will not pursue this problem further in this paper. All remaining values 
of Xdf fr° m data used in the fit are listed in Table 0] Each of these lie around or below 
unity. Since an excellent fit is obtained to DELPHI, SLD and TPC heavy quark-tagged data, 
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we conclude that our fitted heavy quark FFs are reliable even though using them in Eq. 
(fTSJl leads to a poor description of the OPAL heavy quark tagging probabilities. Since the 
DELPHI, SLD and TPC light quark-tagged data is well fitted with the light quark tagging 
probabilities, Eq. is sufficient for describing the latter data. The values of Xbf f° r the 
data to be used for comparison, which were discussed at the beginning of Section llll[ are 
also shown. The serious disagreement with the ALEPH HQ and OPAL Q data found 
here was also found in Ref . |3| , where it was argued that this data has a sizeable contribution 
from charged particles other than the three lightest charged hadrons. For the ALEPH data 
without flavour separation, this argument is supported by the fact that the data for charged 
hadron production significantly overshoots the sum of the hadron identified data. 

In Figs. H3 El we show all these normalized differential cross section data used for fitting 
and for comparison, together with the corresponding theoretical curves from the fit. The 



TPC flavour separated data |17|, particularly the uds quark-tagged data, lie far from their 
theoretical predictions. However, it must be understood that these data are rather old 
compared to the rest of the data used in the fit. At any rate, using them has not affected 
the overall quality of the fit since their errors are large, which explains why their Xbf varues 
in Tableware not too far from unity. Qualitatively, at least, the rise in the calculated cross 
section at low x for decreasing is confirmed by the TPC data, as was first noted in Ref. 



4] 



These figures show that the only TPC data which can significantly constrain a s (Mz) 
are the 7r =l= and identified data shown in Fig. EJ In Fig. we show the gluon-tagged 
three-jet data together with the theoretical curves for D g (x, 4£'j 2 et ). The resulting Xbf values 
shown in the last two lines of Table |U are very high, but it must be kept in mind that the 
theoretical calculation is only correct at LO, and the gluon is only determined at LO. In 
Ref. fsj, where this data was used in the fit, this identification was made only because the 
gluon FF is much less constrained by the remaining data than the quark FFs. 

In Fig. |H1 we compare the longitudinal cross section with the data without flavour sepa- 



from DELPH 
given in Ref. 



32]. The x space coefficient functions of the longitudinal cross section are 
However, since our cross sections are calculated in Mellin space, we cal- 
culate the Mellin transform of these quantities as detailed in Appendix [X] (An alternative 
procedure would be to evolve the FFs in Mellin space as before, and perform the convolution 
of the coefficient functions with the evolved FFs in x space. However, this procedure is nu- 
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merically very slow.) In the unpolarized cross sections used in our fit, the gluon FF for each 
hadron enters only at NLO and so is only determined to LO in our analysis, while it enters 
at LO in the longitudinal cross section, for which a gluon FF determined to NLO is therefore 
required. Thus the curves in Fig. |H] are not completely NLO, but serve to determine the 
quality of our gluon FF. The agreement is excellent for the ALEPH and OPAL data, and 
good for the DELPHI data. Our curves are also very similar to those obtained in Ref. jjj, 
where the LO curves from these authors' LO analysis are also shown. These latter curves 
do not agree with the ALEPH and OPAL data as well as the NLO ones. Thus treating the 
LO gluon FF obtained from their and our fits as NLO results in no loss of consistency in 
this case. 

Finally, we compare cross sections calculated using our FFs for particle production in 
proton-(anti)proton initiated processes with experimental data. Such processes are highly 
dependent on the individual light quark flavour FFs, due to the partonic structure of the 
proton. We use the coefficient functions for the processes a + b — > c + X, where a, b and c 
denote partons, to NLO as calculated in Ref. 35]- We convolute these with our evolved FFs 



for parton c, and the evolved CTEQ6M parton distribution functions [3J| for a and b. Since 

(5) (5) 

our fitted result of A™ =221 MeV is very similar to the result of ALL = 226 MeV obtained 

— MS J MS 



in Ref. 36j, we use the former result in the calculation of a s (fi 2 ). We take Mj = kpj<- The 
cross section at xt = Iprjyfs depends on the FFs for the whole region Xt < z < 1. Since 
we do not (reliably) determine the FFs below z — 0.1 and/or Mf = Mo, we take them in 
this region to be equal to their values at this point. Graphically, we found no discernible 
difference between the resulting predictions and those obtained when the FFs in this region 
were fixed to zero. Firstly, we calculate the invariant differential cross section for inclusive 
7T° production for the process p + p —>■ tt° + X as measured by PHENIX at y/s = 200 GeV 



in Ref. 



□ 



. For this we assume the relation 



Df(x,M}) = ±Df(x,M}) (21) 



to be true, which follows from SU(2) flavour symmetry for pions (see Ref. |37|). Here, is 
the average of the FFs for the processes a, a — > it . (Recall D 1 ^ (x, Mj) is also averaged over 
a and a, but summed over 7r + and ir~ .) The results are shown in Fig. El for k =1/4, 1 and 4, 
together with the PHENIX data. In addition, we also compare the cross section calculated 
from the FFs obtained in Ref. [sj]. For px > 7 GeV, the curve for k — 1 lies closer to the 
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centre of the data than the KKP curve does. Secondly, we calculate the invariant differential 
cross section for inclusive K s production for the process p + p — > K s + X as preliminarily 
measured by STAR at yfs = 200 GeV Ql , and for the process p + p — > K s + X as 
measured by UAl at = 630 GeV in Ref. (40 ]. For this we assume the relation 

D«hx,M 2 f )= 1 -D« ± (x 1 Mj) (22) 

to be true, where b — u, d if a — d, u, otherwise b = a. Eq. (J22)) follows from SU(2) 
flavour symmetry for kaons (see Ref. |2Z(), and is confirmed by the fact that the OPAL 

mesons agree within their errors. 
The predictions are shown in Fig. in a format similar to Fig. El For px > 1.5 GeV, the 
k = 1 curve agrees better with the STAR data than the KKP curve. This disagreement in 
the latter case was observed in Ref. 



|. However, for the older UAl data our predictions 
differ considerably over the whole range, although they are consistent with the data within 
the theoretical errors for p? > 4.5 GeV, while the KKP curve gives good agreement. Since 
the most important difference between our analysis and that of Ref. jb] is the inclusion of 
the OPAL tagging data in our fit, we conclude that this agreement of the KKP curve is 
accidental. 



V. CONCLUSIONS 

This work is an update of the KKP analysis 3], the main difference being that the OPAL 
results on light quark tagging probabilities have been used to phenomenologically constrain 
the individual light quark FFs for the first time. We find that the inclusion of this data 
in the fit makes an important difference to the description of the d,s — *> K ± transitions. 
Light flavour separated FFs are essential for making predictions for inclusive cross sections 
in which there is at least one proton in the initial state and one light hadron in the final 
state (or more than one, in which case other non perturbative quantities are also required 
for subprocesses in which multiple hadrons are emitted from a single parton). Such cross 
sections will be measured, for example, at the LHC In addition, we have included the 



flavour separated TPC data 17[ at y/s = 29 GeV, but such data makes little difference to 
the fit. We have excluded all charged data to be confident that none of the data sets used 
were contaminated with charged particles other than the three lightest charged hadrons. 
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However, good agreement with much of the available charged hadron data, in particular 
that from DELPHI and SLD, was achieved. We point out that although our gluon FF for 
each hadron has been formally determined to LO only, treating it as NLO leads to good 
agreement with the measured longitudinal cross sections in the literature. Finally, relative 
to the KKP predictions, we obtain with our FFs a shift towards the PHENIX data for the 
invariant differential cross section for inclusive ir° production and towards the STAR data 
for the invariant differential cross section for inclusive K s production. 

A determination of a s (Mz) has been performed. We have also calculated the theoretical 
error and find it to be negligible relative to the experimental error. We obtain a s (Mz) = 
0.1176^Q;QQgg, which agrees with the Particle Data Group's world average of a s (M z ) = 
0.1187 ±0.002 H- 

In order to make predictions, our fitted FFs over the range 0.1 < z < 1 
and M < M f < 200 GeV can be obtained from the FORTRAN routines at 
http://www.desy.de/~simon/AKK2005FF.html, which are calculated using cubic spline in- 
terpolation on a linear grid in (z,\nMj). 

APPENDIX A: APPENDIX 

In this appendix, we give all information needed to calculate the longitudinal coefficient 
functions to NLO in Mellin space. 
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The coefficient functions for the longitudinal cross section are given to NLO by 



Cl,ns(z, a a ) = a s C F + a 2 



C 2 F { 451,2(1 - z) - 12Li 3 (-z) + 41nzLi 2 (-z) 



+ 4C(2) ln(l + z) + 4 [21n(l + z)Li 2 (-z) + lnzln 2 (l + z) + 2Sx, 2 {-z)] 
- 4C(2) ln(l - z) - 2 In 2 z ln(l + z) - 4((3) + Q + z\ ln(l - z) 

C(2)" 



-4+ — z~ 2 -8z- -z 3 
5 5 



-3Li 2 (l-z) 



Id : In! ] - : ) ( 4 - 8z - j|z 3 J C(2) + ln 2 (l - z) + f-| + 4z + ^z 3 ) In 2 z 



17 12 



G 



10 5 Z ~ 1 + i Z+ 5 z2)lnZ 



147 12 



— z + -z' 
10 5 



+ C A C F ^ - 2Si, 2 (l - z) - 2 [2 ln(l + z)Li 2 (-z) + lnz ln 2 (l + z) + 25i >2 (-z)] 
+ 6Li 3 (-z) - 2((2) ln(l + z) + 2((2) ln(l - z) - 2 In zLi 2 (-z) + In 2 z ln(l + z) 



„ 6 

2 - -z" 
5 



4 

• 4z + -z' 
5 



+ C(2) ( 4z + -z 3 



2z 



- : ! ) In 2 z - f ln(l - z) + (-g + |^ + §, - ^z 2 ) lnz + 2 C (3) 



1729 6 , 49 4 , 
-z - — z z J 



180 



30 



'2 25 2 

+ C F T R n } \ - (ln(l - z) +lnz) - — + -z 



( 28 

C L ,ps(z,a s ) = a 2 s C F T R n f l 4Li 2 (l - z) +41nzln(l - z) + 6 In 2 z - y - 4z _1 

+ h l z - Az 2 + Qz- 1 - 4z + ^z 2 ) ln(l - z) + (-8 + yz" 1 - 8z + ^z 2 ^) lnz j 
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and 



Ch :g (z, a s ) = ci s Cf 
2 



2 

- - 2 

z 



CU + 2Li 2 (l - z) + 21nzln(l - z) + —C,^)^ 

15 



:■! - ^z 3 ) In 2 z + [ -~ + ^z~ 2 + ^z 3 



+ (-3 + 4z~ 1 -2z) ln(l- 

3 12 31 4 

+ 5"T^ + iT~15 



4 8 2 4 
-3 + 5* +15 
1 12 



i 28 4 2 

2 lnz 



5 5 Z 15 Z 15 



l in 1,2 C(2) 

-9(z)+-bx z __ 



(A3) 



- -Li 2 (l-2) 



+ (-2 + 2Z" 1 ) ln 2 (l - z) - (6 + 8Z" 1 ) In 2 z + ^18 - yV 1 + 2z - j^ 2 J ln(l - z) 
+ C(2) (-8 + 12Z" 1 ) + |l4- M.-1+4Z- V) lnz -81nzln(l - z) 



40 56 , 20 4 2 

1 z 1 z + -z 2 

3 3 3 3 



where the polylogarithms Li n for n — 2, 3, the harmonic sum Si 2 and the function $ are 
defined as 



U 2 (x) 


=_ r ^ 


Lis (a:) 




Jo y 


5i i2 (x) 




2./o y 




rl/(l+x) 




Jx/(l+x) 



(A4) 



In 



To calculate the Mellin transform of the coefficient functions we require only the results in 
Table 1X1 which are obtained from Ref. 4l|. Formally, r\ = (— l) ra , although to analytically 



continue the results in the right hand columns to complex n requires taking 

T] = exp [— 7rlm(n)] [cos (7rRe(n)) + i sin (7rRe(n))] . (A5) 

The harmonic sums Sj(n) are defined for integer n by 

n 1 



Sj{n) 



k=i 



(A6) 
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For complex n, the harmonic sums with j = 1, 2, 3 can be calculated using the results |4 

Si(n) =if)(n + 1) + Je, 

S 2 (n) =-^(n + l) + C(2) and (A7) 
SiCn) =i^(n+l)+C(3), 

where ip(n) and its derivatives can be evaluated for large n using 

1 1 1 1 1 1 / 1 \ 

Hn) _ ln n _____ + __ _ __ + __ _ __ + U\-^j. (A8) 

As noted in Ref. 43^ the harmonic sums when n is small can be calculated by using Eq. 
(IA6|) to write Sj(n) in the form 

Sj {n) = S j (n + r)-j2jj r ^> ( A9 ) 

where r is chosen such that Re(n + r) is large enough to calculate Sj(n + r) using Eq. (jA8jl . 
The function S'(n) (also known as S'_2,i(w)) is defined for integer n by 

fe=l r=l 

A method for calculating S(n) is given in Ref. and generalized for other functions 
appearing in perturbative QCD up to two loop order in Ref. ^_]- Here we present an 
alternative method which is similar to the calculation of the Sj(n) using Eqs. (jA8|) and 
(|A9jl . Firstly, we analytically continue Eq. (jAlOjl to complex values of n by writing it in the 
form 



= E ^ E - - E + «)■ ( An ) 

k=l r=l k=l V ' 

The first term gives — |C(3). For all values of n except for n = —1, —2, the second term 
converges, but very slowly. Instead, we use Eq. (|A10J) to write S(n) in the form 

S(n) = S(n + r) - (-1)" g Jj^^Sfr + k), (A12) 

where r is chosen such that Re(n + r) is large, and calculate S(n + r) as a series in l/(n + r). 
For this purpose, we write Si(n + k) in the form 

oo 1 I 

S 1 {n + k) = lnn + lE + J2-_] Al ™ km > ( A13 ) 



n 

=1 m=0 
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where the A m may be easily calculated using the first relation in Eq. (|A7|) and Eq. (|A8 



then we expand the second term in Eq. (jAll|) in making use of the relation 

d \ ' / I 



£(-l)** r = (x 

k=l \ 



dx 



1 + x 



1 



x=l 



to obtain the result 

, p, k+n 



- (_l)fc ^1 



r=l 



= 0nn + , E) ±J- 2{ - m p + D (4J (j^-l) 

p=l r=0 m=0 v 7 v 

The coefficient of the terms l/n p may now be evaluated, and we find 



x=l 



x=l 



S(n) = -~C(3) 



T](n] 
1 



(Inn + j E ) 



5 11 

+ tt; 7 + 



1 1 

+ 
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2n 2 2n 3 2n 5 2n 7 2n 9 
151 



O 4r 



11 



469 331 67379 „ / 1 

+ T7^+ n +Q "777 



2n 3 12n 4 24n 5 240n 6 240n 7 126n 8 5040n 9 



n 



10 



(A14) 



(A15) 



(A16) 



Since all occurrences of rf must be replaced by unity in the analytic continuation, in the 
second row of Table lA"| we have made the necessary adjustments to that result presented in 
Ref. 0. 

FORTRAN routines for the longitudinal coefficient functions to NLO in Mellin space are 
provided at http://www.desy.de/~simon/cf_long.html. 
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TABLE 1: Values and errors of N, a and j3 in Eq. (|17|) resulting from the 



Hadron 


Flavour 


N 


a 







d 


0.833 ± 0.012 


-1.17 ± 0.01 


1.39 ± 0.02 




u 


0.447 ± 0.007 


-1.58 ± 0.01 


1.01 ± 0.02 




s 


0.519 ± 0.035 


-0.365 ± 0.066 


1.96 ± 0.10 




c 


1.56 ± 0.03 


-1.03 ± 0.01 


3.58 ± 0.07 




b 


0.139 ± 0.001 


-2.24 ± 0.01 


2.77 ± 0.05 




9 


429 ± 3 


2.00 ± 0.01 


5.82 ± 0.01 


K± 


d 


2245 ± 465 


4.14 ± 0.18 


12.0 ± 0.5 




u 


10.9 ± 0.7 


1.72 ± 0.08 


3.44 ± 0.08 




s 


0.529 ± 0.012 


-0.787 ± 0.027 0.915 ± 0.027 




c 


2.28 ± 0.09 


-0.488 ± 0.028 


3.79 ± 0.09 




b 


1.13 ± 0.03 


-0.960 ± 0.016 


6.22 ± 0.09 




9 


15.9 ± 0.5 


2.72 ± 0.05 


2.45 ± 0.03 


p/p 


d 


146 ± 22 


2.30 ± 0.12 


10.4 ± 0.4 




u 


0.0182 ± 0.0014 


-2.37 ± 0.05 


0.507 ± 0.125 




s 


1859 ± 648 


6.67 ± 0.49 


9.17 ± 0.53 




c 


12.0 ± 1.4 


0.860 ± 0.089 


7.50 ± 0.28 




b 


1571 ± 103 


2.19 ± 0.04 


19.0 ± 0.3 




9 


0.867 ± 0.023 


1.13 ± 0.06 


0.854 ± 0.020 
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TABLE 2: Xpw values obtained from the measured light quark tagging probabilities rfc at yfs = 91.2 
GeVinRef. A. 



\. h 

a \. 


7T± 




p/p 


d 


5.05 


0.47 


2.16 


u 


4.87 


0.43 


1.20 


s 


2.69 


0.92 


1.23 



TABLE 3: As in Tabled but for the heavy quarks. 



\^ h 

a 


7T± 




p/p 


b 


17.9 


10.9 


7.64 


c 


24.1 


11.8 


2.96 
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TABLE 4: CM energies, types of data, and Xdf varues f° r various data samples. Samples not used 
in the fits are marked by asterisks, ({h} refers to a sum over light charged hadrons and {q} refers 
to a sum over all 5 flavours of quarks.) 



[GeV] 


Data type 


Xdf 


29.0 


p{h] 

uds 




3.44 [17]* 
2.56 [12* 
1.74 [17]* 






F W} 


0.80 [16] 






uds 


1.01 [17] 








2.51 [17J 








2.14 [17] 






F K 

{<?} 


0.37 [16] 




F p 

{d\ 


0.80 [16] 


91.2 


p{h} 


2.61 [14]* 105 [19]* 
4.99 [8]* 


22.0 [20]* 




jpW 
f r j i 
{uds} 


1.10 [JJ* 64.8 [!£]* 


2.08 [2Q]* 




F {H} 


34.4 [18]* 


0.57 [20]* 






0.21 [14]* 183.6 [18]* 


5.90 [20]* 




pn 
(n\ 


0.98 ri5] 1.13 [14] 


1.82 [8] 




{uds} 


1.82 [14] 


1.12 [8] 




^ C 




1.08 [8] 




F£ 


0.40 [14] 


0.67 [8] 






0.52 [15] 0.31 [14] 


0.52 [8] 




r {uds} 


0.31 [14] 


0.83 [8J 




Ff 




1.79 [8] 




F b K 


0.10 [14] 


1.17 [81 




F p 


0.65 [15] 0.11 [14] 


0.69 [8] 




F p 

{uds} 


0.22 [14] 


1.43 [8] 




PP 
x c 


26 


0.74 [8] 






0.52 PI 41 


1 .14 [81 



Tic 



10" 



10" 



10"' 




10" 



10" 



10"' 



10 



10" 



10"' 



10" 




K 

Tld 



10" 
10" 1 
10 2 

io- 3 

10 4 



— AKK 

KKP 

Kretzer 




K 



10" 
10" 
10 2 

io" 



10 



10" 



10"' 



10"' 



0.5 





10" 

10 

10~ ; 
10' 

10" 
10 ; 
10~ : 
10' 

p 

s 

10" 
10~ ; 
10~ : 
10 ' 



11 






FIG. 1: Light q uarh probabilities fc, .) at * = 91.2 GeV. The dashed curves are calculated 
using the FFs obtained in Ref. |3|, the dotted curves are calculated from the (x,Mj) grid of FFs 
obtained from the analysis of Ref. P] (in which no p/p FFs are obtained), and the solid curves 
are calculated using the FFs obtained in the analysis of this paper. The corresponding measured 
OPAL probabilites of Ref. are also shown. 
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FIG. 3: Normalized differential cross section of inclusive hadron production. The curves are 
calculated from the FFs obtained in our analysis, at 29 (dashed line) and 91.2 (solid line) GeV. 
The upmost, second, third and lowest curves refer to charged hadrons, ir^, and p/p respectively. 
The differential cross section for the charged hadron curve was calculated by ta king the sum of the 
differential cross sections for the three lightest charged hadrons. The ALEPH 15], DELPHI 
OPAL fl, SLD 1 and TPC data sets are shown. The charged hadron data are shown just 
for comparison, but were not used in the fit. Each curve or pair of curves and the corresponding 
data is rescaled relative to the nearest upper one by a factor of 1/5. 
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IG. 4: As in Fig. but for the light quark tagged cross sections. The ALEPH 15], DELPHI 



3, OPAL Q 



, OPAL [20|, SLD || and TPC 



171 ] data sets are shown. The charged hadron data are shown 



just for comparison, but were not used in the fit. 
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FIG. 5: As in Fig. E3 but for the c quark tagged cross sections. The ALEPH Q, OPAL Q , 
SLD [8| and TPC [17| data sets are shown. The two SLD data points at x = 0.654 are for the pion 
(upper) and proton (lower). The charged hadron data are shown just for comparison, but were not 
used in the fit. 
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FIG. 6: As in Fig^H but for the b quark tagged cross sections. The ALEPH Q, DELPHI 
OPAL 20], SLD jy] and TPC data sets are shown. The charged hadron data are shown just 



for comparison, but were not used in the fit. 
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FIG. 7: Gluon FF for charged-hadron production at Mf = 52.4 and 80.2 GeV. The curves are 
calculated from the FFs obtained in our analysis. The three-jet data from ALEPH j^H, with 
Ejet = 26.2 GeV, and from OPAL [22j, with £ jet = 40.1 GeV, are shown. The OPAL data and its 



corresponding curve are rescaled by a factor of 1/100. 
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FIG. 8: Normalized longitudinal differential cross section of inclusive charged hadron production 
at y/s = 91.2 GeV. The curves are calculated from the FFs obtained in our analysis. The data 
sets shown are from ALEPH Q, OPAL Q and DELPHI Q without flavour separation and 
DELPHI [3^1 for light and b quark tagged cross sections. Each curve is rescaled relative to the 
nearest upper one by a factor of 1/30. 
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FIG. 9: The invariant differential cross section for inclusive ir° production in p + p collisions at 
■\fs = 200 GeV. Data from the PHENIX Collaboration [ll| are shown, without the absolute 9.6% 
normalization error. Compared with this data are the cross sections calculated from the FFs 
obtained in this paper (labelled AKK) and that from the FFs of Ref. t 3] (labelled KKP). The 
upper, central and lower AKK curves are calculated with k = 1/4, 1 and 4 respectively. 
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FIG. 10: As in Fig. but for the invariant differential cross section for inclusive K® p roduction in 
pp collisions at y/s = 200 GeV compared with data from the STAR Collaboration [33], and in p+p 
collisions at y/s = 630 GeV compared with data from the UA1 Collaboration 40] . For clarity, the 
former results have been divided by a factor of 30. 
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TABLE A: Mellin transforms required for the Mellin space conversion of the NLO coefficient 
functions of the longitudinal cross section. 
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TABLE A continued: 
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